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Abstract 

Adequate knowledge of the nature of river flow process is crucial for proper 
planning and management of our water resources and environment. This study 
attempts to detect the salient characteristics of flow dynamics of the Karoon 
River in Iran. Daily discharge series observed over a period of six years (1999- 
2004) is analyzed to examine the chaotic and scaling characteristics of the flow 
dynamics. The presence of chaos is investigated through the correlation di- 
mension and Lyapunov exponent methods, while the Hurst exponent and Rcnyi 
dimension analyses are performed to explore the scaling characteristics. The low 
correlation dimension (2.60 ±0.07) and the positive largest Lyapunov exponent 
(0.014 ± 0.001) suggest the presence of low-dimensional chaos; they also imply 
that the flow dynamics are dominantly governed by three variables and can be 
reliably predicted up to 48 days (i.e. prediction horizon). Results from the 
Hurst exponent and Renyi dimension analyses reveal the multifractal character 
of the flow dynamics, with persistent and anti-persistent behaviors observed at 
different time scales. 



1. Introduction 

River flow arises as a result of interactions between climatic inputs and landscape 
characteristics over a wide range of spatial and temporal scales. Depending upon the 
nature of the climatic inputs (e.g. rainfall, temperature) and that of the landscape 
(e.g. basin area, slope, land use), river flow process can have regular or irregular 
behaviors. The degree of regularity or irregularity, i.e. complexity (loosely speaking), 
of the river flow process serves as an important evaluator of its predictability. Study 
on the dynamic nature and predictability of river flow process has always been a key 
research topic in the field of hydrology and water resources and in related fields, as 
such play vital roles both in undertaking short-term water emergency measures and 
in devising long-term water management strategies. 

Linearity and nonlinearity are two of the fundamental properties of river flow pro- 
cess that ultimately dictate its level of complexity. Although the inherent nonlinear 
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nature of river flow process has been known for long [l|, 0, 0] , much of early hydro- 
logic research (especially during 1960s-1980s) was constrained to linear approaches 
0> H EL 01 i the lack of data and computational power largely contributed to this situ- 
ation. Linear approaches continue to be prevalent in river flow and other hydrologic 
studies. However, advances in computational power and data collection during the 
last two decades or so have also facilitated proposal of nonlinear approaches as viable 
alternatives. The nonlinear approaches that have found widespread applications in 
river flow studies include artificial neural networks, data-based mechanistic models 
[||, and chaos theory Q, among others. The present study is concerned with chaos 
theory. In the nonlinear science literature, the term chaos is normally used to refer to 
situations where complex and random-looking behaviors arise from simple nonlinear 
deterministic systems with sensitive dependence on initial conditions [l^] , and the con- 
verse also applies [ill ]. The three fundamental properties inherent in this definition, (i) 
nonlinear interdependence; (ii) hidden determinism and order; and (iii) sensitivity to 
initial conditions, are highly relevant in river systems and processes. For example: (i) 
components and mechanisms involved in the river system (e.g. rainfall, flow, sediment 
transport) act in a nonlinear manner and are also interdependent; (ii) annual cycle in 
river flow possesses determinism and order; and (iii) particle transport in rivers largely 
depend upon the time (i.e. rainy or dry season) at which the particles were released 
at the source, which themselves are often not known. The first property represents 
the 'general' nature of river flow, whereas the second and third represent its 'deter- 
ministic' and 'stochastic' natures, respectively. In view of the obvious relevance of 
its fundamental concepts to river systems, chaos theory has been gaining considerable 
interest in river flow and related studies 0, EH (and Ref. therein). The outcomes of 
these studies are encouraging, as they reveal that chaos theory offers new avenues to 
study the inherent nonlinear and complex nature of the river flow process. 

The major findings in this direction are related to both chaotic and stochastic 
nature in the flow dynamics: it ranges from a less complex (deterministic) to a more 
complex (stochastic) behavior by varying the scale of aggregation. Thus, apparently 
controversial results may emerge by employing analyses on river flow, because of the 
observed transitions from determinism to stochasticity with increasing time scale [lj . 
For an extensive review on chaos theory applications to river flow (and other hydro- 
logic) processes see Ref. [il HR H3| (and references therein). 

Hence, another important aspect in regards to the complexity of river flow process 
is scale. Climatic inputs and landscape characteristics often vary in space and time 
scales and accordingly influence the river processes at various scales. Our general 
perception is that aggregation in scale averages out the variations and reduces the 
level of complexity. Such a perception, however, may not always stand good, since 
averaging may occur only within a limited range of scales, which is often defined by 
the processes themselves. Further, larger spatial and temporal scales may bring their 
own complexities, such as additional terrain types in space and climatic scenarios 
in time. Our knowledge and experience indicate that, for example, daily flow in a 
small river basin often exhibits a higher level of complexity than that of monthly 
flow, but the opposite is often the case when the basin size is large. The last few 
decades have witnessed numerous studies into the scaling properties of river networks 
and river flows [H, Q (and Ref. therein). It is clear, from the above observations, 
that reliable assessment of the complexity and predictability of river flow process and 
identification of the appropriate models for predictions and engineering applications 
requires careful investigation of both the nonlinear (especially chaotic) and scaling 
properties. However, a close examination of the literature reveals that studies have, in 
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general, investigated only the chaotic property or the scaling property, but not both. 

Linear analysis, by means of autocorrelation and power spectrum, and nonlinear 
analysis in the framework of multifractal theory have been extensively used as standard 
methods to investigate only the scaling feature of river flow and rainfall. It has been 
shown that models known as multiplicative cascades, are able to simulate river flow 
scaling behaviour [!?]]. The range and the nature of scaling have been studied on 
several streamflow data from the real world, identifying basic regimes, scale breaking, 
universal multifractal parameters and their independence on basin size, reflecting the 
space-time multiscaling of both the runoff and the rainfall processes [Isl . . Hence, 
space-time multifractals appear to be the natural framework for analyzing the scaling 
features of geodynamical processes including river flows, but it is worth remarking that 
all of these features have been investigated by using river flow and rainfall data. In fact, 
the nature of multiplicative cascades is controversial. Multifractal cascades have been 
successfully introduced as simple stochastic mechanisms for understanding the self- 
similar and intermittent behaviour of turbulent processes [I(|[2l|]. However, it has been 
shown that the deterministic variant of multiplicative cascade models can be chaotic, 
preserving the scaling feature When, as in our case, there is no information about 
rainfall in the same period of river flow data, a deeper comprehension of the dynamics 
requires the analysis of chaotic features as well as multifractal ones: this provides the 
motivation for the present study to investigate both these properties in river flow. 

In this paper, first, we review the main methods for investigating the chaotic be- 
havior and the scaling properties from the time series of a process. The investigation 
on chaos is made by employing a multi-dimensional phase space reconstruction, using 
the embedding theorem [23|, [24| , to obtain preliminary information on possible pat- 
terns and extent of complexity; the correlation dimension method [25| and the largest 
Lyapunov exponent method [26j] are used to investigate the geometry of the phase 
space. The scaling properties are examined through the Hurst exponent estimation 
EJHsl an d the Renyi dimension analysis [2{J. Second, as a case study, flow dynam- 
ics of the Karoon River in Iran is investigated and daily flow data, observed during 
1999-2004, are analyzed by applying these methods. 



2. Methods 

In this study, the investigation of the presence/absence of chaos and scaling be- 
haviors in the river flow series is made using a host of methods. We chose to use 
more than one method to avoid spurious results related to possible drawbacks of each 
procedure. Phase space reconstruction, correlation dimension, and Lyapunov expo- 
nent methods are employed for detecting chaos, whereas Hurst exponent and Renyi 
dimension analyses are performed to identify scaling characteristics. Brief descriptions 
of these methods are presented in the following. 

2.1. Phase space reconstruction 

Phase space is a useful tool for representing the evolution of a system in time. 
Phase space is essentially a graph or a coordinate diagram, whose coordinates rep- 
resent the variables to completely describe the state of the system at any time [231 ] . 
The trajectories in the phase space describe the evolution of the system from some 
initial state, which is assumed to be known, representing the history of the system. 
The region of attraction of these trajectories in the phase space provides important 
qualitative information on the extent of complexity of the system. The idea behind 
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such a reconstruction is that a (nonlinear) dynamic system is characterized by self- 
interaction, and a time series of a single variable can carry the information about the 
dynamics of the entire multivariable system (When multiple variables representing 
the system are available, it is obviously desirable to reconstruct the phase space using 
such multiple variables). Many methods are available for phase space reconstruction 
from a scalar time series. Among these, the method of delays [24f is the most widely 
used one. According to this method, for a scalar time series {x n }, (n — 1, 2, 3, iV) 
multi-dimensional phase space can be reconstructed as 



where m is the dimension of the phase space, called embedding dimension, and r is 
an appropriate (integer) delay time. 

Many different guidelines have been offered in the literature for the selection of 
m and r. For instance, according to the embedding theorem of Takens [U, an m = 
2Do + 1-dimensional phase space is required to completely characterize a dynamic 
system with an attractor dimension D%, whereas, in practice, m > D2 would be 
sufficient Similarly, for the selection of r, some studies suggest the autocorrelation 
function method [U, while mutual information method [U and correlation integral 
method (3^| are also suggested. As of now, there is no general consensus on the 
selection of m and r. In this study, the average mutual information (AMI) method is 
used to select r and the false nearest neighbor (FNN) algorithm is employed to obtain 
m following Ref. [34|, where the best values for the embedding dimension and the lag 
time have been identified within these approaches among the available ones. 

The average mutual information (AMI) approach gets the optimal delay time r as 
the first local minimum [321 of the information measure 



where pi and Pj are the individual probabilities of x n and x n + T respectively, and pij 
is the joint probability. The AMI quantifies the amount of information about x n + T 
if x n is known: when x n + T carries the maximum information about the knowledge of 
x n , AMI is locally minimum. The optimal embedding dimension m can be obtained 
from the false nearest neighbors (FNN) search in phase space [35| : the number of false 
neighbors in the phase space generally changes between two successive embedding 
dimensions mo and mo + 1; the optimal embedding is realized when this number is 
zero. However, real time series are noisy and, hence, the dimension mo + 1 is an 
optimal embedding when the percentage of false neighbors respect to the embedding 
mo is less than a certain threshold, generally fixed to be 1%. 

2.2. Correlation dimension method 

Grassberger and Procaccia [25| introduced a fractal measure, namely the correla- 
tion sum, to quantify the amount of correlations in a time series. Their function is 
defined as the number of pairs of points closer than a given distance e, respect to some 
norm | ■ |, in the embedding space. An improved definition of the correlation sum, 
excluding n time-correlated neighbors in the phase space, is given by 
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where G is the Heaviside step function, n is called Theiler window and the sum is 
extendend to all pairs of points in the embedding space. Correlation sum C(e, m) 
behaves as e m (for any e) for stochastic systems and as t d (for values of e in the scaling 
region) for deterministic ones and some types of coloured noise 0, [3^, [33] ■ When a 
scaling relationship between C(e, m) and e exists, the scaling exponent, namely, the 
correlation dimension D2(e,m), is defined as 

D 2 (e,m)= dl °f } ^ (4) 
aloge 

For Takens embedding theorem, correlation dimension is an estimation of the degrees 
of freedom of the underlying process, i.e. it quantifies the number of equations needed 
to describe the phenomenon if it is deterministic. Numerically, Di(e,rri) can be ob- 
tained by plotting In C(e, m) versus lne for different m and e, and by taking the local 
slopes: for values of m greater than the optimal embedding needed for phase space 
reconstruction, D2 is expected to reach a plateau in the scaling region. The local 
slopes approach is often subjected to poor estimation of the correlation dimension: if 
a scaling region exists, the Takens-Theiler estimator, defined by 



can be used to improve the accuracy of estimation. Chaotic dynamical systems takes 
a fractional value for D2 while deterministic, but not chaotic, dynamical systems 
takes an integer value. For a wide range of stochastic dynamics, D2 goes to infinity: 
unfortunately the existence of a plateau and of a finite and fractional correlation 
dimension is not enough to distinguish chaotic time series from some stochastic ones, 
as coloured noise [3a, |37|] . 



2.3. Lyapunov exponent method 

Processes can be characterized by their sensitivity to initial conditions. If S(t) is 
the distance between two nearby orbits in the phase space at some time t, the evolution 
of the sensitivity to initial conditions £(t) = S(t)/8(0) shows quite different behaviors 
for deterministic and stochastic dynamical systems [381 ] : 

rf£Det.(t) 

dt 
dt 

where A maa! is called largest Lyapunov exponent and H is called Hurst exponent. Recent 
unifying approaches, strictly connected to nonextensive thermodynamics [3Sj, are still 
under investigations. 

Largest Lyapunov exponent is a robust measure of the mean convergence (or di- 
vergence) rate of £(t): for processes that exhibit chaotic behavior, nearby trajectories 
diverge along time and Xmax > 0, while X ma x < for deterministic, but not chaotic, 
dynamical systems. Thus, a positive largest Lyapunov exponent is a strong signature 
of chaos. 

A positive X ma x defines a prediction horizon t*(a) = N*At, i.e. the maximum 
number of samples N* that can be predicted at a sampling time At within an uncer- 
tainty a. Real time series are affected by a measurement error e: it follows that <5(0) = 



Amax^Det. (£) 
,11 



(6) 
(7) 
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e. The worst uncertainty on the prediction of the series is a = max{i„} — min{x n }, i.e. 
it equals the whole range of the series. A desirable uncertainty should be a — 1.96e, 
corresponding to 95% of confidence band. From eq. (0 it follows 

i» = T^log7 (8) 

and t*(1.96e) = X^ ax log 1.96. The largest Lyapunov exponent can be estimated from 
an observed times series with different approaches 0, Eol Elt ]: in the following we 
will adopt Rosenstein's one In the embedded space, the distance of a reference 

point s no from all the other points s n i in its neighborhood of size e, is calculated and 
averaged over the number of neighbors |W S „ | . This procedure, iterated for each point 
along an orbit of AT samples, defines the stretching factor S(e, m, t) that depends on 
the mean convergence (or divergence) rate of nearby trajectories: 



At l ^ At 
S(e,m,t) = — ^2 lo S 



■Ei 



(9) 



where t = AfAt. The largest Lyapunov exponent is the slope of the linear region 
obtained by plotting S(e, m, t) versus t, by keeping fixed e and m. The above procedure 
avoids the estimation of the tangent map, it is fast and easy to implement and it is 
suitable for small and noisy data sets [26l |. However, the estimation of Lyapunov 
exponents may produce spurious results if long time series of high quality are not 
available [42^ or in presence of a strong stochastic contaminating signal [i^ . 

The estimation of the prediction horizon, as previously introduced, can be verified 
by making use of a popular forecasting technique based on nonlinear prediction [4^.|45|. 
The procedure is as follows. In the reconstructed phase space, find all the embedding 
states s„ in the neighborhood W £ (sjv) of size e of the current state sjv. The future 
states s„ 0+ fe, k steps ahead, of all states s„ £ W e (sjv) are successively used for the 
prediction of the measurement at time N + k: 

SN+k = w^r\ Jr Cs) Sno+ * (10) 

i.e. the forecasting is obtained by averaging over all closest embedding states [46|. In 
the presence of an underlying chaotic dynamics, the forecasting error is expected to 
exponentially increase with the forecasting time at a rate \ max , corresponding to the 
largest Lyapunov exponent. The forecast error, divided by the standard deviation of 
the time series, is maximum when the forecast time equals the prediction horizon. 



2.4- Hurst exponent method 

Hurst exponent, previously defined in eq. (ff)l. has been originally introduced for 
investigating the diffusion features of the Nile river (2^, [2j| , and it is widely used to 
detect the scaling regions and to characterize the persistence of a process. Let {y n } 
be the partial sum time series of the series {x n }, defined as 

n ^ N 

i=l 3=1 
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If v is an integer delay time, the Hurst exponent H (q) is obtained from the structure 
function 



t/At— K 

S H {q,u) = _ - ^ |2M +„-2M o r (11) 

' "0=1 

when Sii{q, v) ~ v qIi ^: the power-law is typical of a fractal process at the time scale 
uAt. If H(q) is not a constant function of q, the process is said to be multifractal [38| . 
Hurst exponent is a bounded measure (0 < H(q) < 1) characterizing the persistence 
of a process: 

• H (q) m 0.5 indicates a memoryless time series, with neither short-tem nor long- 
term correlation between states, typical of uncorrelated stochastic processes as 
white noise; 

• < H(q) < 0.5 indicates anti-persistence: increasing trends will be followed 
by decreasing ones, or viceversa, and this behavior tends to be dominant for 
H ->0; 

• 0.5 < H(q) < 1 indicates persistence: there is only one persistent trend typical 
of processes where diffusion is faster than simple brownian motion. 

For a review about Hurst analysis application to hydrological sciences we refer to 
[47| . However, the above method could not be robust if applied to nonstationary 
signals showing evident linear or seasonal trends. Instead, methods based on detrended 
fluctuation analysis (DFA) have been developed to remove linear (or higher order) 
trends, that may exists in nonstationary signals, before performing the scaling analysis. 
We refer to Ref. |48j | for a study of the impact of nonstationary contaminating signals 
on the scaling features of the original data and to Ref. f\j\ for a multifractal study of 
river flow data by means of DFA. Of course, in the case of stationary signals, Hurst 
analysis and DFA agree on the scaling parameters. 

2.5. Generalized q—th order entropy 

Given a scalar time series {x n } defined on some set T>, let us cover this set with 
a partition Ve of disjoint boxes of size e. Let p ( (x) be the probability distribution 
function of the series: p e {xi) is the probability that the series takes the value Xi for 
the partition V t . 

A measure of the average information, needed to specify a point with accuracy e, 
is the Shannon entropy [4j| 



Ui(V e ) = - (logp e (x)) = -^2p e (xi)logp e (xi) 

i 

When a scaling relationship exists between the amount of information and the accuracy 
e, the scaling exponent, called information dimension, is defined as 

Dl = lim 

£ ^o log I 

In a similar way, the concept of entropy and dimension can be respectively gener- 
alized to the q— order Renyi entropy [29| 

H q (V t ) = ^—logJ2[P^)]' J (12) 
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and the q— th order Renyi dimension 



D q = lim (13) 

e^O log I 

where q G [-co, oo] and Hi and Di are obtained if g > 1. The behavior of D q versus 
q defines the fractal spectrum: the underlying process of the time series {x n } is fractal 
when D q is a constant function of q, otherwise it is said to be multifractal. D q can be 
numerically estimated by plotting W q (Ve) versus log - and by taking the slope of the 
linear region. This type of analysis requires long time series in order to avoid spurious 
results. 



3. Data analysis and results 

The Karoon river, with a watershed area of 58, 180 km 2 , is located in southwest of 
the I.R. of Iran, the Khuzestan province is chosen for this study. The river lies between 
the city of Ahwaz (31°20' N,48°4l' E) and the Bahmanshir river (30°25' N,48°12' E), 
which is about 190 km long. The Karoon river is a meandering river which supplies 
water for the irrigation of sugarcane cultivation projects, as well as other agricultural 
lands. 

River flow data, observed over a period of 6 years (from January 1999 to December 
2004), are considered. Fig. [T] shows the variations of daily river flow time series for a 
sampling time of 1 day. 




o I 1 1 1 1 
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Figure 1: Daily Karoon River flow from January 1999 to December 2004. 



Optimal parameters for the phase space reconstruction, namely delay time r = 115 
days and embedding dimension m = 9, are obtained from the first local minimum of 
AMI and FNN search, respectively. Correlation dimension D-^im, e) is estimated from 
the Takens-Theiler measure ([5} up to an embedding dimension m = 18: Fig. [2] shows 
a plateau for D2 = 2.60±0.07, a necessary, but no sufficient, condition for the evidence 
of low-dimensional chaotic dynamics, as previously discussed. 

The stretching factor (j9j), needed for the estimation of the largest Lyapunov 
exponent, is shown in Fig. [3] We found a positive exponent, namely X ma x = 
(0.014 ± 0.001) day -1 , and a prediction horizon of 48.1 ± 3.4 days for a = 1.96e. 

We verified the estimation of the prediction horizon by making use of the technique 
based on forecasting, previously described. In Fig. [3] (inner panel) is shown the 
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Figure 2: Takcns-Thcilcr estimator for the correlation dimension D2(m,e) and plateau fit 
(inner panel) corresponding to D2 = 2.60 ± 0.07 for m > 7 (dashed line). 
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Figure 3: Stretching factor and largest Lyapunov exponent, \ m ax = (0.014 ± 0.001) day -1 , 
corresponding to the slope of the dashed line. Inner panel. Normalized forecast error versus 
the forecast time: the maximum error is obtained for t = 46 days; the prediction horizon 
obtained from the Lyapunov analysis corresponds to 48.1 ± 3.4 days. 



normalized forecast error versus the forecast time: the maximum error is obtained for 
t = 46 days, corresponding to a prediction horizon in excellent agreement with that 
one estimated from Lyapunov analysis. It is worth remarking that our estimation of 
both largest Lyapunov exponent and prediction horizon, characterize global features 
of the underlying dynamics, because they are evaluated by averages on hundreds of 
reference embedding states. 

The structure function 1)111) is estimated for several values of q at some different 
time scales, as shown in Fig. [4] (upper panel): scaling, and thus multifractal behavior, 
emerges at those scales, according to a recent analysis with a multifractal detrended 
fluctuation method on different river flow data [B(J. Hurst exponents, for each q, are 
estimated from the slopes of the straight lines in Fig. [4] (upper panel): a significant 
slope change defines a time scale braking. In Fig. [4](lower panel) Hurst exponents H(q) 
are shown versus q. We identified three significant time scales from slope variations: 

1. From 1 to 28 days; 

2. From 29 to 60 days; 

3. From 61 to 114 days. 




At time scales 1) and 2), Hurst exponents depend on q, as for multifractal dynamical 
systems, H (q) > 0.5 for all q £ [0.5, 4.5] and thus the time series shows a persis- 
tent behavior typical of processes exhibiting long-range dependence. At time scale 3) 
multifractality is stronger than previous time scales, and a transition emerges from 
persistent to anti-persistent behavior through a memoryless state around q = 3, as 
shown in Fig. [4] (lower panel). 

Finally, the q— th order Renyi dimension D q is estimated for several values of 
q £ [0.5,4.5], as shown in Fig. [5] multifractal behavior emerges again, according to 
the Hurst exponent analysis, from an information theoretic point of view. 

Recent studies [Ell . [53 . [53 ] reveal that complex behaviors as multifractality, self- 
organized criticality and on-off intermittency can emerge from nonlinearly-filtered lin- 
ear autoregressive processes (NFLAP). For a correct interpretations of our results, we 
follow the procedure suggested in Ref. [Hi], |H2, 53]. First, we generate 1000 surrogates 
with the same length of the Karoon river flow time series; second, for each surrogate 
time series, we estimate the correlation dimension, the largest Lyapunov exponent, 
the Renyi dimensions and the Hurst exponents. Two types of surrogates, for testing 
two different null hypothesis, are obtained as follows: 

• First type: a surrogate is generated by shuffling the data through a random per- 
mutation of {x n }. The shuffled time series has the same probability distribution 
function of the original one, but time correlations are completely destroyed; 

• Second type: a surrogate is generated with the same probability distibution func- 
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Figure 5: The q— th order Renyi dimension D q estimated from generalized Renyi entropy 
(inner panel) for 0.5 < q < 4.5. 



tion and (almost) the same power spectrum of {x n }. The amplitude-adjusted 
Fourier transform (AAFT) algorithm [H^, [HJ is used to obtain this surrogate. 

First type surrogates are used to test the data against the null hypothesis of an 
underlying uncorrelated stochastic process; second type surrogates are used to test 
the data against the null hypothesis of an underlying linear process distorted by a 
nonlinear filter. 

We found no surrogates, of first or second type, showing, at the same time, a low 
fractional correlation dimension, a positive largest Lyapunov exponent and multifractal 
behavior from both Renyi and Hurst analyses. Thus, we can reject both the null 
hypotheses of an underlying uncorrelated stochastic process or an underlying linear 
process distorted by a nonlinear filter. 

4. Discussion 

In the last years a strong interest emerged for the underlying dynamics of river 
flow. Many studies reveal that both chaotic and stochastic behavior may emerge 
at different spatial and temporal scales. Indeed, multifractality appears to be an 
important feature of this process. However, each study focused or on chaotic features 
either on multifractal ones, but not both at the same time. In fact, deterministic and 
stochastic multiplicative cascades, mainly adopted as models to explain data, show 
similar multifractal signatures and it is not possible to deduce the real nature of the 
river flow without making use of both chaos and scaling analyses. 

Within the present work, we investigated the salient characteristics of dynamics of 
the Karoon river (Iran), by examining the daily discharge time series over a period of 
six years (1999-2004). We followed a nonlinear approach to detect the chaotic and the 
scaling characteristics of the flow dynamics: the presence of chaos has been analyzed 
through the correlation dimension and largest Lyapunov exponent methods, while the 
scaling features have been explored through the Hurst and Renyi analyses. 

Both fractional correlation dimension (2.60 ± 0.07) and positive largest Lyapunov 
exponent (0.014±0.001) suggest the presence of low-dimensional chaos: flow dynamics 
are dominantly governed by three degrees of freedom and can be reliably predicted 
up to 48 days. The estimation of the prediction horizon is in excellent agreement 
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with the value obtained from the nonlinear forecasting analysis: the forecast error 
increases with the forecast time and it reaches its maximum for a time scale of 46 
days. According to recent studies, our results reveal the presence of scaling typical of 
(chaotic) deterministic dynamical systems, although the apparently irregular behavior 
of the data. The fractal structure of a strange attractor emerges in the phase space 
and, of consequence, the underlying dynamics of the Karoon river can be successfully 
modelled by few deterministic equations. However, in the absence of a realistic model 
of the river flow, future discharge can be only monthly predicted from past and current 
measurements. Because of our lack of information on rainfall and atmospherical data 
in the same region, we can not directly quantify their effect on the prediction horizon, 
although our results are in good agreement with recent studies employing different 
forecasting techniques [H^, [571 . IBH , [Ha . loot . 

Results from the Hurst exponent analysis avoid a memoryless phenomenon and re- 
veal, at different time scales, persistent or anti-persistent behavior of the river flow. In- 
deed, nor a unique Hurst exponent neither a unique Renyi dimension can be attributed 
to the entire process, suggesting that river discharge is characterized by anomalous 
scaling typical of multifractal dynamics. In particular, Hurst analysis puts in evidence 
three time scaling regimes. The first and the second ones, corresponding to monthly 
and bimonthly scales, show a persistent behavior: long-range dependence dominates 
the underlying dynamics and diffusion is faster than a simple brownian motion. The 
third scale, ranging from 60 to about 115 days, is the most extended one: for q < 3 
the process is still persistent, whereas the transition from persistent to anti-persistent 
behavior is evident for q > 3. Anti-persistence is symptomatic of a diffusion pro- 
cess slower than a standard brownian motion: in the case of Karoon river flow, it 
emerges in last scaling regime depending on the way we look at the data by means 
of q. Interestingly, 115 days corresponds to the same time that minimizes the aver- 
age mutual information, quantifying the maximum temporal delay, between different 
measurements, before both can be considered no more correlated from an information 
theoretic point of view. Unfortunately, because of our lack of further data, we are not 
able to directly relate these results to rainfall or seasonal effects. The dependence on 
q of the Renyi generalized dimensions is another important signature of multifractal 
behavior, although it does not produce useful information on the real nature of the 
underlying dynamics. Our findings from scaling analyses agree, in general, with recent 
results on the investigation of the multifractal nature of the river flow. 

Finally, we performed the same analyses on two types of surrogate time series, to 
test the null hypotheses that Karoon river flow is an uncorrelated stochastic process 
or a linear stochastic process distorted by a nonlinear filter. We found no surrogates 
showing, at the same time, similar chaotic and scaling characteristics of the flow dy- 
namics and we thus rejected both the null hypotheses. 
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